1 CONFIGURATIONS

# Set the working directory and tool paths on your local computer.
WD <- '/Users/Alec/motrpac/20210826_pass1c-umich'
# Set the gsutil path

knitr::opts_knit$set(root.dir=WD)
knitr::opts_chunk$set(echo = TRUE)
knitr::opts_chunk$set(warning = FALSE)
knitr::opts_chunk$set(message = FALSE)
knitr::opts_chunk$set(cache = FALSE)
1.0.0.0.1 CSS Top Styling
write_css = FALSE
if(write_css){
  writeLines("td, th { padding : 3px } th { background-color:white; color:black; border:1px solid black; text-align:center } td {color:black; border:1px solid black; word-wrap:break-word; white-space:nowrap; overflow: hidden; text-overflow: ellipsis; max-width:300px; text-align:left}", con=file.path(normalizePath(WD), "style/style.css"))
}
1.0.0.0.2 Overview
1.0.0.0.3 Goals of Analysis
  • Examine data dimensions
  • Examine zero, negative, and missing Values
  • Impute missing values
  • Visualize sample-by-sample correlations
  • Identify outliers
  • Identify major causes of variance (drift, sex, control group, hour)
  • Visualize sample-by-feature heatmaps under different transformations
  • Examine sample median and sd distributions between each transformation
  • Transform data
  • Remove outlier samples
  • Save the processed data

2 Prepare Environment & Load Data

2.1 Setup the Environment

################################################################################
##### Resources and Dependencies ###############################################
################################################################################
# Whether to knit document and display data
knit_time = TRUE

# Load dependencies
pacs...man <- c("tidyverse","kableExtra","devtools","MotrpacBicQC","impute","glue",
                "rethinking")
for(pac in pacs...man){
  suppressWarnings(suppressPackageStartupMessages(library(pac, character.only = TRUE)))
}

#browseVignettes("MotrpacBicQC")
############################################################
##### Functions ############################################
############################################################

# Name functions
select <- dplyr::select
map <- purrr::map
desc <- dplyr::desc
arrange <- dplyr::arrange
melt <- reshape2::melt
mutate <- dplyr::mutate
glue <- glue::glue

# Global options
options(dplyr.print_max = 100)
options(scipen=10000)

# Colors
redblue100<-rgb(read.table(paste0(WD,'/colors/redblue100.txt'),sep='\t',row.names=1,header=T))

2.2 Log Variables (1)

ds <- 'pass1a'
site <- 'umich'
tech <- 'rpneg'
TIS <- list('PLA', 'HIP', 'GAS', 'HRT', 'KID', 'LUN', 'LIV', 'BADI', 'WADI')
tis <- list('pla', 'hip', 'gas', 'hrt', 'kid', 'lun', 'liv', 'badi', 'wadi')

2.3 Load Phenotype data

# Phenotype Data (1A)
#########################
pheno_file <- glue("{WD}/data/20201021_pass1a-06-pheno-viallabel_steep.txt")
pheno_df1a <- read.csv(pheno_file, header = T, sep = '\t')
pheno_df1a <- pheno_df1a %>%
  mutate(tissue = case_when(Specimen.Processing.sampletypedescription == 'Brown Adipose' ~ 'badi',
                            Specimen.Processing.sampletypedescription == 'EDTA Plasma' ~ 'pla',
                            Specimen.Processing.sampletypedescription == 'Gastrocnemius' ~ 'gas',
                            Specimen.Processing.sampletypedescription == 'Heart' ~ 'hrt',
                            Specimen.Processing.sampletypedescription == 'Kidney' ~ 'kid',
                            Specimen.Processing.sampletypedescription == 'Liver' ~ 'liv',
                            Specimen.Processing.sampletypedescription == 'Lung' ~ 'lun',
                            Specimen.Processing.sampletypedescription == 'White Adipose' ~ 'wadi',
                            Specimen.Processing.sampletypedescription == 'PaxGene RNA' ~ 'pax',
                            Specimen.Processing.sampletypedescription == 'Hippocampus' ~ 'hip',
                            Specimen.Processing.sampletypedescription == 'Cortex' ~ 'cor',
                            Specimen.Processing.sampletypedescription == 'Hypothalamus' ~ 'hyp',
                            Specimen.Processing.sampletypedescription == 'Vastus Lateralis' ~ 'vas',
                            Specimen.Processing.sampletypedescription == 'Tibia' ~ 'tib',
                            Specimen.Processing.sampletypedescription == 'Aorta' ~ 'aor',
                            Specimen.Processing.sampletypedescription == 'Small Intestine' ~ 'sma',
                            Specimen.Processing.sampletypedescription == 'Adrenals' ~ 'adr',
                            Specimen.Processing.sampletypedescription == 'Colon' ~ 'col',
                            Specimen.Processing.sampletypedescription == 'Spleen' ~ 'spl',
                            Specimen.Processing.sampletypedescription == 'Testes' ~ 'tes',
                            Specimen.Processing.sampletypedescription == 'Ovaries' ~ 'ova'))
pheno_df1a$viallabel <- as.character(pheno_df1a$viallabel)

2.4 Load Data

2.4.1 Load the sample order files

# #created in 20210910_pass1a-umich-sample-annotation_steep.Rmd
# order_file <- glue("{WD}/data/20200910_pass1a1c-sample-order_steep.txt")
# sample.order<-read.delim(order_file,header=T, sep="\t")

# Load the prior pass1a data (takes a few minutes)
pass1a_nested_file <- glue("{WD}/../20200915_metabolomics-pass1a/data/20201010_pass1a-metabolomics-countdata-nested_steep.rds")
pass1a_df <- readRDS(pass1a_nested_file)

# pla
pla_rpneg_pass1a.0.order <- pass1a_df %>%
  filter(STUDY_INSTITUTE == 'University of Michigan') %>%
  filter(NAMED == 'named') %>%
  filter(DATASET == 'untargeted') %>%
  filter(METAB_FAMILY == 'rpneg') %>%
  filter(TISSUE == 'plasma') %>%
  select(SAMPLE_DATA) %>%
  unnest(SAMPLE_DATA) %>%
  left_join(y = pheno_df1a, by = c('sample_id' = 'viallabel')) %>%
  select(sample_id, sample_type, sample_order, Key.anirandgroup, Registration.sex)

# hip
hip_rpneg_pass1a.0.order <- pass1a_df %>%
  filter(STUDY_INSTITUTE == 'University of Michigan') %>%
  filter(NAMED == 'named') %>%
  filter(DATASET == 'untargeted') %>%
  filter(METAB_FAMILY == 'rpneg') %>%
  filter(TISSUE == 'hippocampus') %>%
  select(SAMPLE_DATA) %>%
  unnest(SAMPLE_DATA) %>%
  left_join(y = pheno_df1a, by = c('sample_id' = 'viallabel')) %>%
  select(sample_id, sample_type, sample_order, Key.anirandgroup, Registration.sex)

# gas
gas_rpneg_pass1a.0.order <- pass1a_df %>%
  filter(STUDY_INSTITUTE == 'University of Michigan') %>%
  filter(NAMED == 'named') %>%
  filter(DATASET == 'untargeted') %>%
  filter(METAB_FAMILY == 'rpneg') %>%
  filter(TISSUE == 'gastrocnemius') %>%
  select(SAMPLE_DATA) %>%
  unnest(SAMPLE_DATA) %>%
  left_join(y = pheno_df1a, by = c('sample_id' = 'viallabel')) %>%
  select(sample_id, sample_type, sample_order, Key.anirandgroup, Registration.sex)

# hrt
hrt_rpneg_pass1a.0.order <- pass1a_df %>%
  filter(STUDY_INSTITUTE == 'University of Michigan') %>%
  filter(NAMED == 'named') %>%
  filter(DATASET == 'untargeted') %>%
  filter(METAB_FAMILY == 'rpneg') %>%
  filter(TISSUE == 'heart') %>%
  select(SAMPLE_DATA) %>%
  unnest(SAMPLE_DATA) %>%
  left_join(y = pheno_df1a, by = c('sample_id' = 'viallabel')) %>%
  select(sample_id, sample_type, sample_order, Key.anirandgroup, Registration.sex)

# kid
kid_rpneg_pass1a.0.order <- pass1a_df %>%
  filter(STUDY_INSTITUTE == 'University of Michigan') %>%
  filter(NAMED == 'named') %>%
  filter(DATASET == 'untargeted') %>%
  filter(METAB_FAMILY == 'rpneg') %>%
  filter(TISSUE == 'kidney') %>%
  select(SAMPLE_DATA) %>%
  unnest(SAMPLE_DATA) %>%
  left_join(y = pheno_df1a, by = c('sample_id' = 'viallabel')) %>%
  select(sample_id, sample_type, sample_order, Key.anirandgroup, Registration.sex)

# lun
lun_rpneg_pass1a.0.order <- pass1a_df %>%
  filter(STUDY_INSTITUTE == 'University of Michigan') %>%
  filter(NAMED == 'named') %>%
  filter(DATASET == 'untargeted') %>%
  filter(METAB_FAMILY == 'rpneg') %>%
  filter(TISSUE == 'lung') %>%
  select(SAMPLE_DATA) %>%
  unnest(SAMPLE_DATA) %>%
  left_join(y = pheno_df1a, by = c('sample_id' = 'viallabel')) %>%
  select(sample_id, sample_type, sample_order, Key.anirandgroup, Registration.sex)

# liv
liv_rpneg_pass1a.0.order <- pass1a_df %>%
  filter(STUDY_INSTITUTE == 'University of Michigan') %>%
  filter(NAMED == 'named') %>%
  filter(DATASET == 'untargeted') %>%
  filter(METAB_FAMILY == 'rpneg') %>%
  filter(TISSUE == 'liver') %>%
  select(SAMPLE_DATA) %>%
  unnest(SAMPLE_DATA) %>%
  left_join(y = pheno_df1a, by = c('sample_id' = 'viallabel')) %>%
  select(sample_id, sample_type, sample_order, Key.anirandgroup, Registration.sex)

# badi
badi_rpneg_pass1a.0.order <- pass1a_df %>%
  filter(STUDY_INSTITUTE == 'University of Michigan') %>%
  filter(NAMED == 'named') %>%
  filter(DATASET == 'untargeted') %>%
  filter(METAB_FAMILY == 'rpneg') %>%
  filter(TISSUE == 'brown-adipose') %>%
  select(SAMPLE_DATA) %>%
  unnest(SAMPLE_DATA) %>%
  left_join(y = pheno_df1a, by = c('sample_id' = 'viallabel')) %>%
  select(sample_id, sample_type, sample_order, Key.anirandgroup, Registration.sex)

# wadi
wadi_rpneg_pass1a.0.order <- pass1a_df %>%
  filter(STUDY_INSTITUTE == 'University of Michigan') %>%
  filter(NAMED == 'named') %>%
  filter(DATASET == 'untargeted') %>%
  filter(METAB_FAMILY == 'rpneg') %>%
  filter(TISSUE == 'white-adipose') %>%
  select(SAMPLE_DATA) %>%
  unnest(SAMPLE_DATA) %>%
  left_join(y = pheno_df1a, by = c('sample_id' = 'viallabel')) %>%
  select(sample_id, sample_type, sample_order, Key.anirandgroup, Registration.sex)

rm(pass1a_df)

# Create a list of annotation matrixes
rpneg_pass1a.0.order <- list(pla_rpneg_pass1a.0.order, hip_rpneg_pass1a.0.order, gas_rpneg_pass1a.0.order, 
                 hrt_rpneg_pass1a.0.order, kid_rpneg_pass1a.0.order, lun_rpneg_pass1a.0.order, 
                 liv_rpneg_pass1a.0.order, badi_rpneg_pass1a.0.order, wadi_rpneg_pass1a.0.order)

# Save the Sample order file
save(pla_rpneg_pass1a.0.order, hip_rpneg_pass1a.0.order, gas_rpneg_pass1a.0.order, 
                 hrt_rpneg_pass1a.0.order, kid_rpneg_pass1a.0.order, lun_rpneg_pass1a.0.order, 
                 liv_rpneg_pass1a.0.order, badi_rpneg_pass1a.0.order, wadi_rpneg_pass1a.0.order,
     file = glue("{WD}/data/UM-rpneg/rpneg_pass1a.0.order.RData"))

2.4.2 Load the matrices as .RData files

# UM-rpneg
load(file = glue("{WD}/data/UM-rpneg/UM_rpneg.0.RData"))
pla1a.0 <- pla_rpneg_pass1a.0
hip1a.0 <- hip_rpneg_pass1a.0
gas1a.0 <- gas_rpneg_pass1a.0
hrt1a.0 <- hrt_rpneg_pass1a.0
kid1a.0 <- kid_rpneg_pass1a.0
lun1a.0 <- lun_rpneg_pass1a.0
liv1a.0 <- liv_rpneg_pass1a.0
badi1a.0 <- badi_rpneg_pass1a.0
wadi1a.0 <- wadi_rpneg_pass1a.0

# Create a list of matrices
pass1a.0 <- list(pla1a.0, hip1a.0, gas1a.0, hrt1a.0, kid1a.0, lun1a.0, liv1a.0, badi1a.0, wadi1a.0)

2.5 Remove Internal Standards

is <- list()
for(i in 1:length(pass1a.0)){
  # Remove internal standards
  is[[i]] <- colnames(pass1a.0[[i]])[grepl("istd", colnames(pass1a.0[[i]]), ignore.case = TRUE)]
  # Subset matrix
  pass1a.0[[i]] <- pass1a.0[[i]][, !colnames(pass1a.0[[i]]) %in% is]
}

2.5.1 Create Annotations

i <- 1
tmp.iter1a <- tmp.ref1a <- tmp.sample1a <- pass1a.0.nr <- list()
for(i in 1:length(rpneg_pass1a.0.order)){
  print(tis[[i]])
  # Collect the distances from reference samples
  tmp.iter1a[[i]] <- rpneg_pass1a.0.order[[i]] %>%
    mutate(reference = ifelse(str_sub(sample_id,1,1) == '9', 0, 1)) %>%
    mutate(drift = ifelse(sample_type == 'QC-DriftCorrection', 1, 0)) %>%
    arrange(sample_order)
  tmp.iter1a[[i]]$right_p <- 0
  for(j in 1:nrow(tmp.iter1a[[i]])){
    # set t
    if(tmp.iter1a[[i]][j,'reference'] == 1){
      t = 0
    }else if(tmp.iter1a[[i]][j,'reference'] == 0){
      t = t + 1
    }
    tmp.iter1a[[i]][j,"right_p"] <- t
  }
  t=0
  tmp.iter1a[[i]]$right_p_d <- 0
  for(j in 1:nrow(tmp.iter1a[[i]])){
    # set t
    if(tmp.iter1a[[i]][j,'drift'] == 1){
      t = 0
    }else if(tmp.iter1a[[i]][j,'drift'] == 0){
      t = t + 1
    }
    tmp.iter1a[[i]][j,"right_p_d"] <- t
  }
  t=0

  tmp.iter1a[[i]] <- tmp.iter1a[[i]] %>%
    arrange(desc(sample_order))
  tmp.iter1a[[i]]$left_p <- 0
  for(j in 1:nrow(tmp.iter1a[[i]])){
    # set t
    if(tmp.iter1a[[i]][j,'reference'] == 1){
      t = 0
    }else if(tmp.iter1a[[i]][j,'reference'] == 0){
      t = t + 1
    }
    tmp.iter1a[[i]][j,"left_p"] <- t
  }
  t=0
  tmp.iter1a[[i]]$left_p_d <- 0
  for(j in 1:nrow(tmp.iter1a[[i]])){
    # set t
    if(tmp.iter1a[[i]][j,'drift'] == 1){
      t = 0
    }else if(tmp.iter1a[[i]][j,'drift'] == 0){
      t = t + 1
    }
    tmp.iter1a[[i]][j,"left_p_d"] <- t
  }
  t=0

  tmp.iter1a[[i]] <- tmp.iter1a[[i]] %>%
    rowwise() %>%
    mutate(min_p = min(c(left_p,right_p) ,na.rm = TRUE)) %>%
    mutate(sum_p = sum(c(left_p,right_p) ,na.rm = TRUE)) %>%
    mutate(min_p_d = min(c(left_p_d,right_p_d) ,na.rm = TRUE)) %>%
    mutate(sum_p_d = sum(c(left_p_d,right_p_d) ,na.rm = TRUE)) %>%
    arrange(sample_order)
  tmp.join1a <- tmp.iter1a[[i]] %>%
    select(sample_id, sample_order, left_p, right_p, min_p, sum_p, left_p_d, right_p_d, min_p_d, sum_p_d)

  # Collect the sample order for test+ref samples
  tmp.ref1a[[i]] <- rpneg_pass1a.0.order[[i]] %>%
    arrange(sample_order) %>% 
    mutate(control = ifelse(grepl('Control', Key.anirandgroup), 1, 0)) %>%
    mutate(drift = ifelse(grepl('Drift', sample_type), 1, 0)) %>%
    mutate(reference = ifelse(str_sub(sample_id,1,1) == '9', 0, 1)) %>%
    mutate(time = case_when(grepl('IPE', Key.anirandgroup) ~ 0,
                            grepl('0 hr', Key.anirandgroup) ~ 0,
                            grepl('0.5 hr', Key.anirandgroup) ~ 0.5,
                            grepl('1 hr', Key.anirandgroup) ~ 1,
                            grepl('4 hr', Key.anirandgroup) ~ 4,
                            grepl('7 hr', Key.anirandgroup) ~ 7,
                            grepl('24 hr', Key.anirandgroup) ~ 24,
                            grepl('48 hr', Key.anirandgroup) ~ 48)) %>%
    left_join(y = tmp.join1a) %>%
    select(sample_id,sample_order, Registration.sex, control, time, drift, reference,
         left_p, right_p, min_p, sum_p, left_p_d, right_p_d, min_p_d, sum_p_d)
  N2 <- tmp.ref1a[[i]][,1] %>% unlist()
  miss1 <- N2[!N2 %in% row.names(pass1a.0[[i]])] # Verify all samples are in the pass1a.0[[i]] file
  miss1
  # sample.order %>% filter(sample_id == "90750016606")
  print('Vice Versa:')
  miss2 <- row.names(pass1a.0[[i]])[!row.names(pass1a.0[[i]]) %in% N2]
  miss2
  N2 <- N2[!N2 %in% c(miss1,miss2)] # TODO: investigate why samples are missing, continue for now
  # Reorder pass1a.0[[i]] by run order
  pass1a.0[[i]] <- pass1a.0[[i]][N2,]
  all(N2 == row.names(pass1a.0[[i]])) # Must be true
  tmp.ref1a[[i]] <- tmp.ref1a[[i]] %>%
    filter(sample_id %in% N2) %>%
    select(sample_order, Registration.sex, control, time, drift, reference,
           left_p, right_p, min_p, sum_p, left_p_d, right_p_d, min_p_d, sum_p_d) %>% as.matrix()

  # Collect the sample order for test samples
  options(digits = 14)
  tmp.sample1a[[i]] <- rpneg_pass1a.0.order[[i]] %>%
    filter(substr(sample_id, 1, 1) == '9') %>%
    arrange(sample_order) %>% 
    mutate(control = ifelse(grepl('Control', Key.anirandgroup), 1, 0)) %>%
    mutate(time = case_when(grepl('IPE', Key.anirandgroup) ~ 0,
                            grepl('0 hr', Key.anirandgroup) ~ 0,
                            grepl('0.5 hr', Key.anirandgroup) ~ 0.5,
                            grepl('1 hr', Key.anirandgroup) ~ 1,
                            grepl('4 hr', Key.anirandgroup) ~ 4,
                            grepl('7 hr', Key.anirandgroup) ~ 7,
                            grepl('24 hr', Key.anirandgroup) ~ 24,
                            grepl('48 hr', Key.anirandgroup) ~ 48)) %>%
    left_join(y = tmp.join1a) %>%
    mutate(sample_id = str_replace_all(sample_id, pattern = '-', replacement = '')) %>%
    mutate(sample_id = as.numeric(sample_id)) %>%
    select(sample_id, sample_order, Registration.sex, control, time, 
           left_p, right_p, min_p, sum_p, left_p_d, right_p_d, min_p_d, sum_p_d) %>%
    as.matrix()

  N1 <- row.names(pass1a.0[[i]])[substr(row.names(pass1a.0[[i]]), 1, 1) == '9']
  tmp.sample1a[[i]] <- tmp.sample1a[[i]][tmp.sample1a[[i]][,1] %in% N1,]
  # Reorder pass1a.0[[i]].nr by run order
  pass1a.0.nr[[i]] <- pass1a.0[[i]][N1,]
  tmp.sample1a[[i]][,1][!tmp.sample1a[[i]][,1] %in% row.names(pass1a.0.nr[[i]])] # Verify all samples are in the pass1a.0[[i]] file
  all(as.character(tmp.sample1a[[i]][,1]) == row.names(pass1a.0.nr[[i]]))
  # If out of order, command below will ensure pass1a.0[[i]].nr in run order
  #pass1a.0[[i]].nr <- pass1a.0[[i]].nr[as.character(tmp.sample1a[[i]][,1]),]
}
## [1] "pla"
## [1] "Vice Versa:"
## [1] "hip"
## [1] "Vice Versa:"
## [1] "gas"
## [1] "Vice Versa:"
## [1] "hrt"
## [1] "Vice Versa:"
## [1] "kid"
## [1] "Vice Versa:"
## [1] "lun"
## [1] "Vice Versa:"
## [1] "liv"
## [1] "Vice Versa:"
## [1] "badi"
## [1] "Vice Versa:"
## [1] "wadi"
## [1] "Vice Versa:"

3 Dimensions, Zero/Neg/Missing Values, & Log2

3.1 Dimensions (with reference samples)

# Lists
NR <- P <- list()
for(i in 1:length(pass1a.0)){
  print(tis[[i]])
  NR[[i]] <- dim(pass1a.0[[i]])[1]
  P[[i]] <- dim(pass1a.0[[i]])[2]
  print(dim(pass1a.0[[i]]))
}
## [1] "pla"
## [1] 96 97
## [1] "hip"
## [1] 113 139
## [1] "gas"
## [1] 98 78
## [1] "hrt"
## [1] 120  94
## [1] "kid"
## [1] 120  87
## [1] "lun"
## [1] 120 102
## [1] "liv"
## [1] 103  81
## [1] "badi"
## [1] 116 143
## [1] "wadi"
## [1] 112  78

3.2 Dimensions (without reference samples)

N <- list()
for(i in 1:length(pass1a.0.nr)){
  print(tis[[i]])
  N[[i]] <- dim(pass1a.0.nr[[i]])[1]
  print(dim(pass1a.0.nr[[i]]))
}
## [1] "pla"
## [1] 77 97
## [1] "hip"
## [1]  78 139
## [1] "gas"
## [1] 78 78
## [1] "hrt"
## [1] 78 94
## [1] "kid"
## [1] 78 87
## [1] "lun"
## [1]  78 102
## [1] "liv"
## [1] 78 81
## [1] "badi"
## [1]  78 143
## [1] "wadi"
## [1] 76 78

3.3 Negative or Zero Values

confirmed: no negative or zero values

for(i in 1:length(pass1a.0.nr)){
  print(tis[[i]])
  print(min(pass1a.0[[i]],na.rm=T))
}
## [1] "pla"
## [1] 36.40427817
## [1] "hip"
## [1] 71.37304298
## [1] "gas"
## [1] 49.21625509
## [1] "hrt"
## [1] 120.973489
## [1] "kid"
## [1] 101.3145201
## [1] "lun"
## [1] 61.68128751
## [1] "liv"
## [1] 502.8654934
## [1] "badi"
## [1] 100.201935
## [1] "wadi"
## [1] 63.25290985

3.4 Missing Features (with references)

Blank reference samples at the beginning and end

pass1a.0.f.c0 <- list()
for(i in 1:length(pass1a.0)){
  pass1a.0.f.c0[[i]] <-apply(pass1a.0[[i]],1,function(x) sum(is.na(x))) 
}

# Plot
par(mfrow=c(3,3), mar = c(2,3,2,0))
for(i in 1:length(pass1a.0)){
  plot(pass1a.0.f.c0[[i]], ylim = c(0,max(unlist(P))), main = glue("{tis[[i]]}"), xlab = 'Samples', ylab = 'Missing Values')
}

3.5 Blank Samples (without references)

No blank test samples

pass1a.0.nr.f.c0 <- list()
for(i in 1:length(pass1a.0.nr)){
  pass1a.0.nr.f.c0[[i]] <-apply(pass1a.0.nr[[i]],1,function(x) sum(is.na(x))) 
}

# Plot
par(mfrow=c(3,3), mar = c(2,3,2,0))
for(i in 1:length(pass1a.0.nr)){
  plot(pass1a.0.nr.f.c0[[i]], ylim = c(0,max(unlist(P))), main = glue("{tis[[i]]}"), xlab = 'Samples', ylab = 'Missing Values')
}

3.6 Examine Distribution of missing values (with reference samples)

pass1a.0.f.c0 <- list()
for(i in 1:length(pass1a.0.nr)){
  pass1a.0.f.c0[[i]] <- apply(pass1a.0[[i]],2,function(x) sum(is.na(x))) 
}
# Plot
par(mfrow=c(3,3), mar = c(2,3,2,0))
for(i in 1:length(pass1a.0.nr)){
  plot(pass1a.0.f.c0[[i]], ylim = c(0,max(unlist(NR))), main = glue("{tis[[i]]}"), xlab = 'Features', ylab = 'Missing Values')
}

3.7 Examine Distribution of missing values (test samples only)

pass1a.0.nr.f.c0 <- list()
for(i in 1:length(pass1a.0.nr)){
  pass1a.0.nr.f.c0[[i]] <- apply(pass1a.0.nr[[i]],2,function(x) sum(is.na(x))) 
}
# Plot
par(mfrow=c(3,3), mar = c(2,3,2,0))
for(i in 1:length(pass1a.0.nr)){
  plot(pass1a.0.nr.f.c0[[i]], ylim = c(0,max(unlist(N))), main = glue("{tis[[i]]}"), xlab = 'Features', ylab = 'Missing Values'); abline(h = 20)
}

3.8 Remove high-missing features

Remove high missing features

rm_n <- list()
for(i in 1:length(pass1a.0.nr.f.c0)){
  print(tis[[i]])
  rm_n[[i]] <- sum(pass1a.0.nr.f.c0[[i]]>=20)
  pass1a.0[[i]] <- pass1a.0[[i]][,pass1a.0.nr.f.c0[[i]]<20]
  pass1a.0.nr[[i]] <- pass1a.0.nr[[i]][,pass1a.0.nr.f.c0[[i]]<20]
  print(dim(pass1a.0[[i]]))
  print(dim(pass1a.0.nr[[i]]))
}
## [1] "pla"
## [1] 96 96
## [1] 77 96
## [1] "hip"
## [1] 113 109
## [1]  78 109
## [1] "gas"
## [1] 98 76
## [1] 78 76
## [1] "hrt"
## [1] 120  94
## [1] 78 94
## [1] "kid"
## [1] 120  87
## [1] 78 87
## [1] "lun"
## [1] 120 102
## [1]  78 102
## [1] "liv"
## [1] 103  81
## [1] 78 81
## [1] "badi"
## [1] 116 123
## [1]  78 123
## [1] "wadi"
## [1] 112  76
## [1] 76 76

3.9 Take the log2

pass1a.0.1 <- pass1a.0.nr1 <- list()
for(i in 1:length(pass1a.0)){
  pass1a.0.1[[i]] <- log(pass1a.0[[i]], 2)
  pass1a.0.nr1[[i]] <- log(pass1a.0.nr[[i]], 2)
}

3.10 Total missing values (includes reference samples)

for(i in 1:length(pass1a.0.1)){
  print(glue("{tis[[i]]}"))
  print(sum(is.na(pass1a.0.1[[i]])))
}
## pla
## [1] 2
## hip
## [1] 408
## gas
## [1] 20
## hrt
## [1] 268
## kid
## [1] 289
## lun
## [1] 363
## liv
## [1] 105
## badi
## [1] 563
## wadi
## [1] 450

3.11 Total missing values (test samples)

feature_impute <- list()
for(i in 1:length(pass1a.0.nr1)){
  print(tis[[i]])
  print(sum(is.na(pass1a.0.nr1[[i]])))
  feature_impute[[i]] <- apply(is.na(pass1a.0.nr1[[i]]),2,sum)[apply(is.na(pass1a.0.nr1[[i]]),2,sum) > 0]
}
## [1] "pla"
## [1] 0
## [1] "hip"
## [1] 3
## [1] "gas"
## [1] 20
## [1] "hrt"
## [1] 19
## [1] "kid"
## [1] 9
## [1] "lun"
## [1] 5
## [1] "liv"
## [1] 80
## [1] "badi"
## [1] 10
## [1] "wadi"
## [1] 45

3.12 Confirm New Distribution of missing features (test samples only)

pass1a.0.nr1.f.c0 <- list()
for(i in 1:length(pass1a.0.nr1)){
  pass1a.0.nr1.f.c0[[i]] <- apply(pass1a.0.nr1[[i]],2,function(x) sum(is.na(x))) 
}

# Plot
par(mfrow=c(3,3), mar = c(2,3,2,0))
for(i in 1:length(pass1a.0.nr1)){
  plot(pass1a.0.nr1.f.c0[[i]], ylim = c(0,max(unlist(NR))), main = glue("{tis[[i]]}"), xlab = 'Features', ylab = 'Missing Values'); abline(h = 20)
}

3.13 Log Variables (2)

# NR
# N
# P
neg_vals <- 0
zero_vals <- 0
feature_na_filter <- 20
knn_k <- 10


P1 <- NR1 <- N1 <- na_vals_impute <- list()
for(i in 1:length(pass1a.0)){
  print(tis[[i]])
  P1[[i]] <- dim(pass1a.0.nr[[i]])[2]
  NR1[[i]] <- dim(pass1a.0[[i]])[1]
  N1[[i]] <- dim(pass1a.0.nr[[i]])[1]
  na_vals_impute[[i]] <- sum(is.na(pass1a.0.nr1[[i]]))
  print(feature_impute[[i]])
}
## [1] "pla"
## named integer(0)
## [1] "hip"
## glutamate-13C5 [iSTD] 
##                     3 
## [1] "gas"
## Tetradecadienoic acid      Myristoleic acid         Deoxycytidine 
##                     3                     2                     1 
##     Pentadecylic acid            Xanthosine 
##                     1                    13 
## [1] "hrt"
##        alanine-13C3 [iSTD] 3-Methyl-2-oxovaleric acid 
##                          2                          3 
##                Ketoleucine                Lauric acid 
##                          2                          8 
##            Lignoceric acid       Ursodeoxycholic acid 
##                          3                          1 
## [1] "kid"
##          Lactic acid          Capric acid Pregnenolone sulfate 
##                    2                    4                    2 
##          Cholic acid 
##                    1 
## [1] "lun"
## Tauro-alpha-muricholic acid 
##                           5 
## [1] "liv"
##                      Tetradecadienoic acid 
##                                          2 
##                           Myristoleic acid 
##                                          2 
##                              Myristic acid 
##                                          2 
##                   gamma-Glutamylisoleucine 
##                                          9 
##              alpha-N-Phenylacetylglutamine 
##                                          1 
##                    Octadecatetraenoic acid 
##                                          2 
##                            Eicosenoic acid 
##                                          1 
##                      Docosapentaenoic acid 
##                                          8 
##                        Docosatrienoic acid 
##                                          2 
##                         Tetracosenoic acid 
##                                          2 
## 5alpha-Androstane-3alpha-ol-17-one sulfate 
##                                         18 
##                       Ursodeoxycholic acid 
##                                          6 
##                       Pregnenolone sulfate 
##                                          8 
##                                Cholic acid 
##                                         17 
## [1] "badi"
##    glutamine-13C5 [iSTD]          Lignoceric acid     Ursodeoxycholic acid 
##                        1                        2                        6 
## Taurohyodeoxycholic acid 
##                        1 
## [1] "wadi"
##         alanine-13C3 [iSTD]                 Ketoleucine 
##                           5                           3 
##        Hydroxyoctanoic acid              Kynurenic acid 
##                           1                           1 
##     citric acid-13C6 [iSTD]       N-Acetylphenylalanine 
##                           2                          14 
##      Hydroxydodecanoic acid   Hydroxytetradecanoic acid 
##                           2                           5 
##                  Xanthosine       Heptadecanedioic acid 
##                           1                           2 
##                Behenic acid Tauro-alpha-muricholic acid 
##                           1                           7 
##            Taurocholic acid 
##                           1

4 Imputation

4.1 Imputation (test samples)

pass1a.0.nr2 <- list()
for(i in 1:length(pass1a.0)){
  if(na_vals_impute[[i]] > 0){
    print(tis[[i]])
    glue("Features & Values to impute:")
    feature_impute[[i]]
    pass1a.0.nr2[[i]] <-impute.knn(pass1a.0.nr1[[i]],k=10)$data
    #view the features before and after imputation
    par(mfrow=c(2,1),bg="black")
    image(as.matrix(pass1a.0.nr1[[i]][,pass1a.0.nr1.f.c0[[i]]>0]),col=redblue100,axes=F)
    image(as.matrix(pass1a.0.nr2[[i]][, pass1a.0.nr1.f.c0[[i]]>0]),col=redblue100,axes=F)
    par(mfrow=c(1,1) ,bg="white")
    glue("Verified no missing values: {sum(is.na(pass1a.0.nr2[[i]]))}")  #verified 0
  }else{
    print('No missing values to impute')
    pass1a.0.nr2[[i]] <- pass1a.0.nr1[[i]]
  }
}
## [1] "No missing values to impute"
## [1] "hip"

## [1] "gas"

## [1] "hrt"

## [1] "kid"

## [1] "lun"

## [1] "liv"

## [1] "badi"

## [1] "wadi"

5 NxN heatmaps

5.1 Plotting NxN heatmaps (+ reference samples)

a <- list(0.82,0.85,0.80,
       0.80,0.80,0.82,
       0.80,0.80,0.80)
b <- 1
par(mfrow=c(3,3), mar = c(0,0,3,0))
i <- 1
for(i in 1:length(pass1a.0)){
  print(tis[[i]])
  x <- tmp.ref1a[[i]][,1]
  sidebar  <- ( (b - a[[i]]) * ((x - min(x))/((max(x) - min(x)))) + a[[i]]) 
  sidebar <- cbind(sidebar,sidebar,sidebar)
  cor.tmp<-cor(t(pass1a.0.1[[i]]),method="spearman",use="pairwise.complete.obs") #includes ref samples
  print(dim(cor.tmp))
  image(cbind(cor.tmp,sidebar),
            col=redblue100,axes=FALSE,zlim=c(a[[i]],1), main=glue("{TIS[[i]]}2-{NR1[[i]]},
                                                                  run-order, z=({a[[i]]},1)"), asp = 1)
}
## [1] "pla"
## [1] 96 96
## [1] "hip"
## [1] 113 113
## [1] "gas"
## [1] 98 98
## [1] "hrt"
## [1] 120 120
## [1] "kid"
## [1] 120 120
## [1] "lun"
## [1] 120 120
## [1] "liv"
## [1] 103 103
## [1] "badi"
## [1] 116 116
## [1] "wadi"
## [1] 112 112

# if(knit_time){
#   rpneg_pass1a.0.order[[i]] %>%
#     arrange(sample_order) %>%
#     select(sample_id, sample_type, sample_order) %>%
#     knitr::kable(format = "html") %>%
#     scroll_box(width = "100%", height = "400px")
# }

5.2 Plotting NxN heatmaps (test samples) and Identification of Outlier Samples

cor.tmp <- o.s <- list()
par(mfrow=c(3,3), mar = c(0,0,3,0))
for(i in 1:length(pass1a.0)){
  x <- tmp.sample1a[[i]][,2]
  sidebar  <- ( (b - a[[i]]) * ((x - min(x))/((max(x) - min(x)))) + a[[i]]) 
  sidebar <- cbind(sidebar,sidebar,sidebar)
  cor.tmp[[i]]<-cor(t(pass1a.0.nr1[[i]]),method="spearman",use="pairwise.complete.obs")
  glue("Verified {N[[i]]} test samples: {dim(cor.tmp)[1]}") #verified, =78
  image(
    cbind(cor.tmp[[i]][order(tmp.sample1a[[i]][,2]),],sidebar),
    col=redblue100,axes=FALSE, zlim=c(a[[i]],1), main=glue("{TIS[[i]]}2-{N1[[i]]}, 
                                                           Run-Order, z=({a[[i]]},1)"), asp = 1)
}

# Outlier sample filter
o.f <- list(0.89, 0.915, 0.89,
            0.89, 0.90, 0.90,
            0.90, 0.89, 0.87)
par(mfrow=c(3,3), mar = c(2,3,2,0))
for(i in 1:length(pass1a.0)){
  plot(apply(cor.tmp[[i]][order(tmp.sample1a[[i]][,2]),order(tmp.sample1a[[i]][,2])],1,median), main = glue("{tis[[i]]}"),
       xlab = 'Samples (Run-Order)', ylab = 'Cor Medians'); abline(h=o.f[[i]], col = 'blue')
}

# Determine which samples are outliers
for(i in 1:length(pass1a.0)){
  if(is.na(o.f[[i]])){
    o.s[[i]] <- NA
  }else{
    o.s[[i]] <- colnames(cor.tmp[[i]])[apply(cor.tmp[[i]],1,median)<o.f[[i]]]
  }
  glue("Outlier Samples: {length(o.s[[i]])}")
}
# Examine outliers
glue("Outlier Samples in {TIS[[1]]}:")
## Outlier Samples in PLA:
rpneg_pass1a.0.order[[1]] %>% filter(sample_id %in% o.s[[1]]) %>%
  knitr::kable(format = "html") %>% scroll_box(width = "100%", height = "100%")
sample_id sample_type sample_order Key.anirandgroup Registration.sex
90010013104 Sample 31 Exercise - IPE 1
glue("Outlier Samples in {TIS[[2]]}:")
## Outlier Samples in HIP:
rpneg_pass1a.0.order[[2]] %>% filter(sample_id %in% o.s[[2]]) %>%
  knitr::kable(format = "html") %>% scroll_box(width = "100%", height = "100%")
sample_id sample_type sample_order Key.anirandgroup Registration.sex
90010015206 Sample 101 Exercise - IPE 1
glue("Outlier Samples in {TIS[[3]]}:")
## Outlier Samples in GAS:
rpneg_pass1a.0.order[[3]] %>% filter(sample_id %in% o.s[[3]]) %>%
  knitr::kable(format = "html") %>% scroll_box(width = "100%", height = "100%")
sample_id sample_type sample_order Key.anirandgroup Registration.sex
90136015506 Sample 25 Exercise - IPE 2
90156015506 Sample 53 Exercise - 1 hr 1
90007015506 Sample 57 Exercise - 0.5 hr 2
90010015506 Sample 63 Exercise - IPE 1
90012015506 Sample 77 Exercise - 1 hr 1
90025015506 Sample 85 Exercise - 0.5 hr 2
glue("Outlier Samples in {TIS[[4]]}:")
## Outlier Samples in HRT:
rpneg_pass1a.0.order[[4]] %>% filter(sample_id %in% o.s[[4]]) %>%
  knitr::kable(format = "html") %>% scroll_box(width = "100%", height = "100%")
sample_id sample_type sample_order Key.anirandgroup Registration.sex
90156015807 Sample 63 Exercise - 1 hr 1
90012015807 Sample 68 Exercise - 1 hr 1
90010015807 Sample 108 Exercise - IPE 1
glue("Outlier Samples in {TIS[[5]]}:")
## Outlier Samples in KID:
rpneg_pass1a.0.order[[5]] %>% filter(sample_id %in% o.s[[5]]) %>%
  knitr::kable(format = "html") %>% scroll_box(width = "100%", height = "100%")
sample_id sample_type sample_order Key.anirandgroup Registration.sex
90012015906 Sample 74 Exercise - 1 hr 1
90010015906 Sample 104 Exercise - IPE 1
glue("Outlier Samples in {TIS[[6]]}:")
## Outlier Samples in LUN:
rpneg_pass1a.0.order[[6]] %>% filter(sample_id %in% o.s[[6]]) %>%
  knitr::kable(format = "html") %>% scroll_box(width = "100%", height = "100%")
sample_id sample_type sample_order Key.anirandgroup Registration.sex
90010016606 Sample 69 Exercise - IPE 1
90012016606 Sample 98 Exercise - 1 hr 1
glue("Outlier Samples in {TIS[[7]]}:")
## Outlier Samples in LIV:
rpneg_pass1a.0.order[[7]] %>% filter(sample_id %in% o.s[[7]]) %>%
  knitr::kable(format = "html") %>% scroll_box(width = "100%", height = "100%")
sample_id sample_type sample_order Key.anirandgroup Registration.sex
90010016805 Sample 47 Exercise - IPE 1
90007016805 Sample 79 Exercise - 0.5 hr 2
glue("Outlier Samples in {TIS[[8]]}:")
## Outlier Samples in BADI:
rpneg_pass1a.0.order[[8]] %>% filter(sample_id %in% o.s[[8]]) %>%
  knitr::kable(format = "html") %>% scroll_box(width = "100%", height = "100%")
sample_id sample_type sample_order Key.anirandgroup Registration.sex
90012016906 Sample 31 Exercise - 1 hr 1
90010016906 Sample 33 Exercise - IPE 1
90156016906 Sample 72 Exercise - 1 hr 1
glue("Outlier Samples in {TIS[[9]]}:")
## Outlier Samples in WADI:
rpneg_pass1a.0.order[[9]] %>% filter(sample_id %in% o.s[[9]]) %>%
  knitr::kable(format = "html") %>% scroll_box(width = "100%", height = "100%")
sample_id sample_type sample_order Key.anirandgroup Registration.sex
90011017005 Sample 35 Exercise - 1 hr 2
90010017005 Sample 72 Exercise - IPE 1
90156017005 Sample 82 Exercise - 1 hr 1
90012017005 Sample 88 Exercise - 1 hr 1

5.3 Visualize the Sex-Time-Control NxN Heatmap

par(mfrow=c(3,3), mar = c(0,0,3,0))
for(i in 1:length(pass1a.0)){
  x <- tmp.sample1a[[i]][,3]
  sex.type  <- ( (b - a[[i]]) * ((x - min(x))/((max(x) - min(x)))) + a[[i]]) 
  x <- tmp.sample1a[[i]][,4]
  control.type  <- ( (b - a[[i]]) * ((x - min(x))/((max(x) - min(x)))) + a[[i]]) 
  x <- as.factor(tmp.sample1a[[i]][,5]) %>% as.numeric()
  time.type  <- ( (b - a[[i]]) * ((x - min(x))/((max(x) - min(x)))) + a[[i]]) 
  sidebar<-cbind(sex.type,sex.type,sex.type, 
               control.type,control.type,control.type, 
               time.type,time.type,time.type)
  image(
    cbind(
    cor.tmp[[i]][order(tmp.sample1a[[i]][,3],tmp.sample1a[[i]][,4],tmp.sample1a[[i]][,5]),
          order(tmp.sample1a[[i]][,3],tmp.sample1a[[i]][,4],tmp.sample1a[[i]][,5])],
   sidebar[order(tmp.sample1a[[i]][,3],tmp.sample1a[[i]][,4],tmp.sample1a[[i]][,5]),]),
    col=redblue100,axes=FALSE, zlim=c(a[[i]],1), main=glue("{TIS[[i]]}2-{N[[i]]}, 
                                                           Sex-Control-Time, z=({a[[i]]},1)"), asp = 1)
}

par(mfrow=c(3,3), mar = c(2,3,2,0))
for(i in 1:length(pass1a.0)){
  plot(apply(cor.tmp[[i]][order(tmp.sample1a[[i]][,3],tmp.sample1a[[i]][,4],tmp.sample1a[[i]][,5]), 
                         order(tmp.sample1a[[i]][,3],tmp.sample1a[[i]][,4],tmp.sample1a[[i]][,5])],1,median), main = glue("{tis[[i]]}"))
}

5.4 NxN Heatmap: Test and Reference Samples by Run Order

cor.tmp1 <- list()
par(mfrow=c(3,3), mar = c(0,0,3,0))
for(i in 1:length(pass1a.0)){
  x <- is.na(tmp.ref1a[[i]][,2]) %>% as.integer()
  sample.type  <- ( (b - a[[i]]) * ((x - min(x))/((max(x) - min(x)))) + a[[i]]) 
  sidebar<-cbind(sample.type,sample.type,sample.type,sample.type,sample.type)
  cor.tmp1[[i]]<-cor(t(pass1a.0.1[[i]]),method="spearman",use="pairwise.complete.obs") #includes ref samples
  image(
    cbind(cor.tmp1[[i]], sidebar),
    col=redblue100,axes=FALSE, zlim=c(a[[i]],1), main=glue("{TIS[[i]]}2, 
                                                         Run-Order (+Ref-type), z=({a[[i]]},1)"), asp = 1)
}

5.5 Log Variables (3)

run_var <- list(0, 0, 0,
                0, 0, 0,
                0, 0, 0)
outlier_sample_n <- list( length(o.s[[1]]),length(o.s[[2]]),length(o.s[[3]]),
                         length(o.s[[4]]),length(o.s[[5]]),length(o.s[[6]]),
                         length(o.s[[7]]),length(o.s[[8]]),length(o.s[[9]]) )
outlier_samples <- o.s
outlier_filter <- o.f

6 N-P Heatmaps, Median and SD Distributions, & Transformations

6.0.1 N-P Heatmaps: Log2, imputed, sample order,

Note: Samples as columns Sidebar: p-values, t-values, and sig (binary) comparing normal test samples and outlier test samples

n.s <- hmo <- list()
par(mfrow=c(3,3), mar = c(0,0,3,0))
for(i in 1:length(pass1a.0)){
  n.s[[i]] <- row.names(pass1a.0.nr2[[i]])[!row.names(pass1a.0.nr2[[i]]) %in% o.s[[i]]]
  hmo[[i]] <- heatmap(pass1a.0.nr2[[i]])$colInd
}

for(i in 1:length(pass1a.0)){
  image(pass1a.0.nr2[[i]][order(tmp.sample1a[[i]][,2]),hmo[[i]]]
      ,col=redblue100,axes=F,main=glue("{TIS[[i]]}2, 
                                       f{P1[[i]]}-HC, Run-Order"), asp = 1)
}

par(mfrow=c(3,3), mar = c(2,3,2,0))
for(i in 1:length(pass1a.0)){
  plot(apply(pass1a.0.nr2[[i]][order(tmp.sample1a[[i]][,2]),hmo[[i]]],1,median), main = glue("{tis[[i]]}"))
}

6.0.2 N-P Heatmap: Log2, imputed, Sex-Control-Time order (2)

Note: pass1a.0.2 and pass1a.0.nr2 represent knn-imputed and log2 Note: Samples as columns

par(mfrow=c(3,3), mar = c(0,0,3,0))
for(i in 1:length(pass1a.0)){
  a[[i]] <- range(pass1a.0.nr2[[i]])[1]
  b[[i]] <- range(pass1a.0.nr2[[i]])[2]
  x <- tmp.sample1a[[i]][,3]
  sex.type  <- ( (b[[i]] - a[[i]]) * ((x - min(x))/((max(x) - min(x)))) + a[[i]]) 
  x <- tmp.sample1a[[i]][,4]
  control.type <- ( (b[[i]] - a[[i]]) * ((x - min(x))/((max(x) - min(x)))) + a[[i]]) 
  x <- as.factor(tmp.sample1a[[i]][,5]) %>% as.numeric()
  time.type  <- ( (b[[i]] - a[[i]]) * ((x - min(x))/((max(x) - min(x)))) + a[[i]]) 
  sidebar<-cbind(sex.type,sex.type,sex.type,sex.type,sex.type, 
               control.type,control.type, control.type, control.type, control.type,
               time.type, time.type, time.type, time.type, time.type)
  image(
    cbind(pass1a.0.nr2[[i]][order(tmp.sample1a[[i]][,3],tmp.sample1a[[i]][,4],tmp.sample1a[[i]][,5]), hmo[[i]] ], 
      sidebar[order(tmp.sample1a[[i]][,3],tmp.sample1a[[i]][,4],tmp.sample1a[[i]][,5]),]),
      col=redblue100,axes=F,main=glue("{TIS[[i]]}2, 
                                    f{P1[[i]]}-HC, Sex-Control-Time"), asp =1)
}

par(mfrow=c(3,3), mar = c(2,3,2,0))
for(i in 1:length(pass1a.0)){
  plot(apply(pass1a.0.nr2[[i]][order(tmp.sample1a[[i]][,3],tmp.sample1a[[i]][,4],tmp.sample1a[[i]][,5]),
     hmo[[i]]],1,median), main = glue("{tis[[i]]}"))
}

6.0.3 Transformations – only for the test samples. This can be revised to remove outlier samples

Strategy is not functionally programed (must revert log transformed back to linear for pass1a.0.3a)

pass1a.0.3a<-pass1a.0.3b<-pass1a.03c<-pass1a.03c2<-pass1a.02b<-pass1a.03d<-pass1a.03d2<-pass1a.02b2<-pass1a.03d3 <- list()
for(i in 1:length(pass1a.0)){
  # pass1a.0.r3a<-pass1a.0.r3b<-pass1a.0.r3c<-pass1a.0.r3c2<-pass1a.0.r2b<-pass1a.0.r3d<-pass1a.0.r3d2<-pass1a.0.r2b2<-pass1a.0.r3d3 <-pass1a.0.2
  pass1a.0.3a[[i]]<-pass1a.0.3b[[i]]<-pass1a.03c[[i]]<-pass1a.03c2[[i]]<-pass1a.02b[[i]]<-pass1a.03d[[i]]<-pass1a.03d2[[i]]<-pass1a.02b2[[i]]<-pass1a.03d3[[i]]<-pass1a.0.nr2[[i]]
  #pass1a.03c3<-pass1a.0.3b2<-pass1a.0.nr2
}

6.0.4 Examine the median and mean distributions (2)

par(mfrow=c(3,3), mar = c(2,3,2,0))
for(i in 1:length(pass1a.0)){
  tmp.s.median <- apply(pass1a.0.nr2[[i]],1, median)
  tmp.s.mean <- apply(pass1a.0.nr2[[i]],1, mean)
  plot(tmp.s.median,tmp.s.mean, asp = 1, main = glue("{tis[[i]]}")); abline(0,1)
}

6.0.5 Examine the median and sd distribution (1)

tmp.f.median <- tmp.f.sd <- list()
par(mfrow=c(3,3), mar = c(2,3,2,0))
for(i in 1:length(pass1a.0)){
  tmp.f.median[[i]] <- apply(2^pass1a.0.nr2[[i]],2, median)
  tmp.f.sd[[i]] <- apply(2^pass1a.0.nr2[[i]],2, sd)
  plot(y = tmp.f.sd[[i]], x = tmp.f.median[[i]],log="xy", main = glue("{tis[[i]]}"))
}

6.0.6 Examine the median and sd distribution (1; w/wo outlier samples)

tmp2.f.median <- list()
par(mfrow=c(3,3), mar = c(2,3,2,0))
for(i in 1:length(pass1a.0)){
  tmp <- 2^pass1a.0.nr2[[i]][n.s[[i]],]
  #dim(tmp)
  tmp2.f.median[[i]] <- apply(tmp,2, median)
  #tmp2.f.sd <- apply(tmp,2, sd)
  plot(tmp.f.median[[i]],tmp2.f.median[[i]],log="xy", main = glue("{tis[[i]]}"))
  #plot(tmp.f.sd,tmp2.f.sd,log="xy")
}

tmp2.f.sd <- list()
par(mfrow=c(3,3), mar = c(2,3,2,0))
for(i in 1:length(pass1a.0)){
  tmp <- 2^pass1a.0.nr2[[i]][n.s[[i]],]
  tmp2.f.sd[[i]] <- apply(tmp,2, sd)
  plot(tmp.f.sd[[i]],tmp2.f.sd[[i]],log="xy",main = glue("{tis[[i]]}"))
}

6.1 Center by median, scale by standard deviation across features (3a)

Verification of the first feature

par(mfrow=c(3,3), mar = c(2,3,2,0))
j <- 148
for(i in 1:length(pass1a.0)){
  for (j in 1:length(tmp.f.median[[i]])) {
    pass1a.0.3a[[i]][,j]<-(2^pass1a.0.nr2[[i]][,j] - tmp.f.median[[i]][j])/ tmp.f.sd[[i]][j]
  }
  plot(2^pass1a.0.nr2[[i]][,1],pass1a.0.3a[[i]][,1], main = glue("{tis[[i]]}")) #spot-check, verified
}

6.1.1 N-P Heatmap: median-centered, sd-scaled imputed, run order (3a)

hmo2 <- list()
for(i in 1:length(pass1a.0)){
  hmo2[[i]] <- heatmap(pass1a.0.3a[[i]])$colInd
}

# Run Order; Original HC
###########################
par(mfrow=c(3,3), mar = c(0,0,3,0))
for(i in 1:length(pass1a.0)){
  image(pass1a.0.3a[[i]][order(tmp.sample1a[[i]][,2]), hmo[[i]]],
    col=redblue100,axes=F,main=glue("{TIS[[i]]}3a, 
                                    f{P1[[i]]}-HC, Run-Order"), asp = 1)
}

# Run Order; New HC
###########################
par(mfrow=c(3,3), mar = c(0,0,3,0))
for(i in 1:length(pass1a.0)){
  image(pass1a.0.3a[[i]][order(tmp.sample1a[[i]][,2]), hmo2[[i]]],
    col=redblue100,axes=F,main=glue("{TIS[[i]]}3a, 
                                    f{P1[[i]]}-HC(3a), Run-Order"), asp = 1)
}

par(mfrow=c(3,3), mar = c(2,3,2,0))
for(i in 1:length(pass1a.0)){
  plot(apply(pass1a.0.3a[[i]][order(tmp.sample1a[[i]][,2]),
     hmo[[i]]],1,median), main = glue("{tis[[i]]}"), ylim = c(-2,2)); abline(h = 0, col = 'blue')
}

6.1.2 Examine the median and sd distribution (2)

tmp.f.median2 <- tmp.f.sd2 <- list()
par(mfrow=c(3,3), mar = c(2,3,2,0))
for(i in 1:length(pass1a.0)){
  tmp.f.median2[[i]] <- apply(pass1a.0.nr2[[i]],2, median)
  tmp.f.sd2[[i]] <- apply(pass1a.0.nr2[[i]],2, sd)
  plot(y = tmp.f.sd2[[i]], x = tmp.f.median2[[i]],log="xy", main = glue("{tis[[i]]}"))
}

6.1.3 Examine the median and sd distribution (2; w/wo outlier samples)

tmp2.f.median2 <- list()
par(mfrow=c(3,3), mar = c(2,3,2,0))
for(i in 1:length(pass1a.0)){
  tmp <-pass1a.0.nr2[[i]][n.s[[i]],]
  #dim(tmp)
  tmp2.f.median2[[i]] <- apply(tmp,2, median)
  plot(tmp.f.median2[[i]],tmp2.f.median2[[i]],log="xy", main = glue("{tis[[i]]}"))
}

tmp2.f.sd2 <- list()
par(mfrow=c(3,3), mar = c(2,3,2,0))
for(i in 1:length(pass1a.0)){
  tmp <- pass1a.0.nr2[[i]][n.s[[i]],]
  tmp2.f.sd2[[i]] <- apply(tmp,2, sd)
  plot(tmp.f.sd2[[i]],tmp2.f.sd2[[i]],log="xy",main = glue("{tis[[i]]}"))
}

6.2 Log2, Center by median, scale by standard deviation across features (3b)

Verification of the first & last feature

for(i in 1:length(pass1a.0)){
  for (j in 1:length(tmp2.f.median[[i]])) {
    pass1a.0.3b[[i]][,j]<-(pass1a.0.nr2[[i]][,j]- tmp.f.median2[[i]][j])/tmp.f.sd2[[i]][j]
  }
}
par(mfrow=c(3,3), mar = c(2,3,2,0))
for(i in 1:length(pass1a.0)){
  plot(pass1a.0.nr2[[i]][,1],pass1a.0.3b[[i]][,1], main = glue("{tis[[i]]}")) #spot-check, verified
}

for(i in 1:length(pass1a.0)){
  plot(pass1a.0.nr2[[i]][,P1[[i]]],pass1a.0.3b[[i]][,P1[[i]]], main = glue("{tis[[i]]}")) #spot-check, verified
}

6.2.1 N-P Heatmap: median-centered, sd-scaled imputed, run order (3b)

hmo3 <- list()
for(i in 1:length(pass1a.0)){
  hmo3[[i]] <- heatmap(pass1a.0.3b[[i]])$colInd
}

# Run Order; Original HC
###########################
par(mfrow=c(3,3), mar = c(0,0,3,0))
for(i in 1:length(pass1a.0)){
  image(pass1a.0.3b[[i]][order(tmp.sample1a[[i]][,2]), hmo[[i]]],
    col=redblue100,axes=F,main=glue("{TIS[[i]]}3b, 
                                    f{P1[[i]]}-HC, Run-Order"), asp = 1)
}

# Run Order; New HC
###########################
par(mfrow=c(3,3), mar = c(0,0,3,0))
for(i in 1:length(pass1a.0)){
  image(pass1a.0.3b[[i]][order(tmp.sample1a[[i]][,2]), hmo3[[i]]],
    col=redblue100,axes=F,main=glue("{TIS[[i]]}3b, 
                                    f{P1[[i]]}-HC(3b), Run-Order"), asp = 1)
}

par(mfrow=c(3,3), mar = c(2,3,2,0))
for(i in 1:length(pass1a.0)){
  plot(apply(pass1a.0.3b[[i]][order(tmp.sample1a[[i]][,2]),
     hmo[[i]]],1,median), main = glue("{tis[[i]]}"), ylim = c(-1.5,1.5)); abline(h = 0, col = 'blue')
}

6.2.2 N-P Heatmap: median-centered, sd-scaled imputed, Sex-Control-Time order (3b)

par(mfrow=c(3,3), mar = c(0,0,3,0))
for(i in 1:length(pass1a.0)){
  a[[i]] <- range(pass1a.0.3b[[i]])[1]
  b[[i]] <- range(pass1a.0.3b[[i]])[2]
  x <- tmp.sample1a[[i]][,3]
  sex.type  <- ( (b[[i]] - a[[i]]) * ((x - min(x))/((max(x) - min(x)))) + a[[i]]) 
  x <- tmp.sample1a[[i]][,4]
  control.type <- ( (b[[i]] - a[[i]]) * ((x - min(x))/((max(x) - min(x)))) + a[[i]]) 
  x <- as.factor(tmp.sample1a[[i]][,5]) %>% as.numeric()
  time.type  <- ( (b[[i]] - a[[i]]) * ((x - min(x))/((max(x) - min(x)))) + a[[i]]) 
  sidebar<-cbind(sex.type,sex.type,sex.type,sex.type,sex.type, 
               control.type,control.type, control.type, control.type, control.type,
               time.type, time.type, time.type, time.type, time.type)
  image(
    cbind(pass1a.0.3b[[i]][order(tmp.sample1a[[i]][,3],tmp.sample1a[[i]][,4],tmp.sample1a[[i]][,5]), hmo[[i]] ], 
      sidebar[order(tmp.sample1a[[i]][,3],tmp.sample1a[[i]][,4],tmp.sample1a[[i]][,5]),]),
      col=redblue100,axes=F,main=glue("{TIS[[i]]}2, 
                                    f{P1[[i]]}-HC, Sex-Control-Time"), asp =1)
}

par(mfrow=c(3,3), mar = c(2,3,2,0))
for(i in 1:length(pass1a.0)){
  plot(apply(pass1a.0.3b[[i]][order(tmp.sample1a[[i]][,3],tmp.sample1a[[i]][,4],tmp.sample1a[[i]][,5]),
     hmo[[i]]],1,median), main = glue("{tis[[i]]}"), ylim = c(-1.5,1.5)); abline(h = 0, col = 'blue')
}

6.2.3 Examine the sample median and sd distributions (3b)

par(mfrow=c(3,3), mar = c(0,2,1,0))
for(i in 1:length(pass1a.0)){
  boxplot(data.frame(t(pass1a.0.3b[[i]])), ylim = c(-2,2), main = glue("{tis[[i]]}"), xaxt = "n")
}

glue("Outlier samples removed:")
## Outlier samples removed:
par(mfrow=c(3,3), mar = c(0,2,1,0))
for(i in 1:length(pass1a.0)){
  boxplot(data.frame(t(pass1a.0.3b[[i]][n.s[[i]],])), ylim = c(-2,2), main = glue("{tis[[i]]}"), xaxt = "n")
}

6.2.4 TODO (in additional script): Collect sample medians independent of run-order outliers to determine if samples should be centered and/or scaled

6.2.5 TODO: Incorporate outlier sample removal

6.2.6 TODO: Incorporate flagged samples

6.2.7 TODO: Incorporate outlier feature removal from transformations

7 Save the processed data

7.1 Concatenate the Processing Decisions

pass1a.0.vars <- list()
comments <- list(NA, NA, NA,
                 NA, NA, NA,
                 NA, NA, NA)
pass1a.0.vars <- data.frame()
for(i in 1:length(pass1a.0)){
  feature_impute[[i]] = names(feature_impute[[i]])
  # The processing decisions
  pass1a.0.vars <- rbind(pass1a.0.vars, data.frame(ds, site, tech, 'tis' = tis[[i]], 'NR' = NR[[i]], 'N' = N[[i]], 'P' = P[[i]], 
                                   neg_vals, zero_vals, feature_na_filter, 'P1' = P1[[i]], 'NR1' = NR1[[i]], 
                                   'N1' = N1[[i]], 'imputed_features' = paste(feature_impute[[i]], collapse = '; '), 
                                   'na_vals_impute' = na_vals_impute[[i]], knn_k, 'run_var' = run_var[[i]], 
                                   'outlier_sample_n' = outlier_sample_n[[i]], 'outlier_samples' =  paste(o.s[[i]], collapse = '; '), 
                                   'outlier_filter' = o.f[[i]],  'internal_standards_n' = length(is[[i]]), 
                                   'internal_standards' = paste(is[[i]], collapse = '; '), 'comments' = comments[[i]]))
}

# visualize the processing decisions
if(knit_time){
  pass1a.0.vars %>%
    knitr::kable(format = "html") %>%
    scroll_box(width = "100%", height = "100%")
}
ds site tech tis NR N P neg_vals zero_vals feature_na_filter P1 NR1 N1 imputed_features na_vals_impute knn_k run_var outlier_sample_n outlier_samples outlier_filter internal_standards_n internal_standards comments
pass1a umich rpneg pla 96 77 97 0 0 20 96 96 77 0 10 0 1 90010013104 0.890 10 valine-13C5 [iSTD]; aspartate-13C4 [iSTD]; glutamine-13C5 [iSTD]; glutamate-13C5 [iSTD]; phenylalanine-13C9 [iSTD]; tyrosine-13C9 [iSTD]; citric acid-13C6 [iSTD]; tryptophan-15N2 [iSTD]; Gibberellic acid [iSTD]; 24-epi-Brassinolide [iSTD] NA
pass1a umich rpneg hip 113 78 139 0 0 20 109 113 78 glutamate-13C5 [iSTD] 3 10 0 1 90010015206 0.915 22 alanine-13C3 [iSTD]; lactate-13C3 [iSTD]; serine-13C3 [iSTD]; valine-13C5 [iSTD]; threonine-13C4 [iSTD]; Asparagine-15N2 [iSTD]; aspartate-13C4 [iSTD]; isoleucine-13C6 [iSTD]; Leucine-13C/Isoleucine-13C6 [iSTD]; leucine-13C6 [iSTD]; malate-13C4 [iSTD]; anthranilic acid-15N [iSTD]; glutamine-13C5 [iSTD]; glutamate-13C5 [iSTD]; lysine-13C6 [iSTD]; histidine-13C6 [iSTD]; phenylalanine-13C9 [iSTD]; tyrosine-13C9 [iSTD]; citric acid-13C6 [iSTD]; tryptophan-15N2 [iSTD]; Gibberellic acid [iSTD]; 24-epi-Brassinolide [iSTD] NA
pass1a umich rpneg gas 98 78 78 0 0 20 76 98 78 Tetradecadienoic acid; Myristoleic acid; Deoxycytidine; Pentadecylic acid; Xanthosine 20 10 0 6 90136015506; 90156015506; 90007015506; 90010015506; 90012015506; 90025015506 0.890 0 NA
pass1a umich rpneg hrt 120 78 94 0 0 20 94 120 78 alanine-13C3 [iSTD]; 3-Methyl-2-oxovaleric acid; Ketoleucine; Lauric acid; Lignoceric acid; Ursodeoxycholic acid 19 10 0 3 90156015807; 90012015807; 90010015807 0.890 12 alanine-13C3 [iSTD]; serine-13C3 [iSTD]; valine-13C5 [iSTD]; aspartate-13C4 [iSTD]; glutamine-13C5 [iSTD]; glutamate-13C5 [iSTD]; phenylalanine-13C9 [iSTD]; tyrosine-13C9 [iSTD]; citric acid-13C6 [iSTD]; tryptophan-15N2 [iSTD]; Gibberellic acid [iSTD]; 24-epi-Brassinolide [iSTD] NA
pass1a umich rpneg kid 120 78 87 0 0 20 87 120 78 Lactic acid; Capric acid; Pregnenolone sulfate; Cholic acid 9 10 0 2 90012015906; 90010015906 0.900 10 alanine-13C3 [iSTD]; valine-13C5 [iSTD]; aspartate-13C4 [iSTD]; glutamine-13C5 [iSTD]; glutamate-13C5 [iSTD]; tyrosine-13C9 [iSTD]; citric acid-13C6 [iSTD]; tryptophan-15N2 [iSTD]; Gibberellic acid [iSTD]; 24-epi-Brassinolide [iSTD] NA
pass1a umich rpneg lun 120 78 102 0 0 20 102 120 78 Tauro-alpha-muricholic acid 5 10 0 2 90010016606; 90012016606 0.900 16 alanine-13C3 [iSTD]; lactate-13C3 [iSTD]; aspartate-13C4 [iSTD]; isoleucine-13C6 [iSTD]; leucine-13C6 [iSTD]; anthranilic acid-15N [iSTD]; glutamine-13C5 [iSTD]; glutamate-13C5 [iSTD]; lysine-13C6 [iSTD]; histidine-13C6 [iSTD]; phenylalanine-13C9 [iSTD]; tyrosine-13C9 [iSTD]; citric acid-13C6 [iSTD]; tryptophan-15N2 [iSTD]; Gibberellic acid [iSTD]; 24-epi-Brassinolide [iSTD] NA
pass1a umich rpneg liv 103 78 81 0 0 20 81 103 78 Tetradecadienoic acid; Myristoleic acid; Myristic acid; gamma-Glutamylisoleucine; alpha-N-Phenylacetylglutamine; Octadecatetraenoic acid; Eicosenoic acid; Docosapentaenoic acid; Docosatrienoic acid; Tetracosenoic acid; 5alpha-Androstane-3alpha-ol-17-one sulfate; Ursodeoxycholic acid; Pregnenolone sulfate; Cholic acid 80 10 0 2 90010016805; 90007016805 0.900 0 NA
pass1a umich rpneg badi 116 78 143 0 0 20 123 116 78 glutamine-13C5 [iSTD]; Lignoceric acid; Ursodeoxycholic acid; Taurohyodeoxycholic acid 10 10 0 3 90012016906; 90010016906; 90156016906 0.890 22 alanine-13C3 [iSTD]; lactate-13C3 [iSTD]; serine-13C3 [iSTD]; valine-13C5 [iSTD]; threonine-13C4 [iSTD]; Asparagine-15N2 [iSTD]; aspartate-13C4 [iSTD]; isoleucine-13C6 [iSTD]; Leucine-13C/Isoleucine-13C6 [iSTD]; leucine-13C6 [iSTD]; malate-13C4 [iSTD]; anthranilic acid-15N [iSTD]; glutamine-13C5 [iSTD]; glutamate-13C5 [iSTD]; lysine-13C6 [iSTD]; histidine-13C6 [iSTD]; phenylalanine-13C9 [iSTD]; tyrosine-13C9 [iSTD]; citric acid-13C6 [iSTD]; tryptophan-15N2 [iSTD]; Gibberellic acid [iSTD]; 24-epi-Brassinolide [iSTD] NA
pass1a umich rpneg wadi 112 76 78 0 0 20 76 112 76 alanine-13C3 [iSTD]; Ketoleucine; Hydroxyoctanoic acid; Kynurenic acid; citric acid-13C6 [iSTD]; N-Acetylphenylalanine; Hydroxydodecanoic acid; Hydroxytetradecanoic acid; Xanthosine; Heptadecanedioic acid; Behenic acid; Tauro-alpha-muricholic acid; Taurocholic acid 45 10 0 4 90011017005; 90010017005; 90156017005; 90012017005 0.870 9 alanine-13C3 [iSTD]; valine-13C5 [iSTD]; aspartate-13C4 [iSTD]; glutamine-13C5 [iSTD]; glutamate-13C5 [iSTD]; citric acid-13C6 [iSTD]; tryptophan-15N2 [iSTD]; Gibberellic acid [iSTD]; 24-epi-Brassinolide [iSTD] NA

7.2 Save the Data

# pla
pla1a.0.nr1 <- pass1a.0.nr1[[1]]
pla1a.0.1 <- pass1a.0.1[[1]]
pla1a.0.nr2 <- pass1a.0.nr2[[1]]
pla1a.0.3a <- pass1a.0.3a[[1]]
pla1a.0.3b <- pass1a.0.3b[[1]]
pla1a.0.vars <- pass1a.0.vars[1,]
# save the data
save(pla1a.0.nr1, pla1a.0.1, pla1a.0.nr2, pla1a.0.3a, pla1a.0.3b, pla1a.0.vars,
  file = glue("{WD}/data/UM-rpneg/UM_rpneg_processed_pla1a.0.RData"))

# hip
hip1a.0.nr1 <- pass1a.0.nr1[[2]]
hip1a.0.1 <- pass1a.0.1[[2]]
hip1a.0.nr2 <- pass1a.0.nr2[[2]]
hip1a.0.3a <- pass1a.0.3a[[2]]
hip1a.0.3b <- pass1a.0.3b[[2]]
hip1a.0.vars <- pass1a.0.vars[2,]
# save the data
save(hip1a.0.nr1, hip1a.0.1, hip1a.0.nr2, hip1a.0.3a, hip1a.0.3b, hip1a.0.vars,
  file = glue("{WD}/data/UM-rpneg/UM_rpneg_processed_hip1a.0.RData"))

# gas
gas1a.0.nr1 <- pass1a.0.nr1[[3]]
gas1a.0.1 <- pass1a.0.1[[3]]
gas1a.0.nr2 <- pass1a.0.nr2[[3]]
gas1a.0.3a <- pass1a.0.3a[[3]]
gas1a.0.3b <- pass1a.0.3b[[3]]
gas1a.0.vars <- pass1a.0.vars[3,]
# save the data
save(gas1a.0.nr1, gas1a.0.1, gas1a.0.nr2, gas1a.0.3a, gas1a.0.3b, gas1a.0.vars,
  file = glue("{WD}/data/UM-rpneg/UM_rpneg_processed_gas1a.0.RData"))

# hrt
hrt1a.0.nr1 <- pass1a.0.nr1[[4]]
hrt1a.0.1 <- pass1a.0.1[[4]]
hrt1a.0.nr2 <- pass1a.0.nr2[[4]]
hrt1a.0.3a <- pass1a.0.3a[[4]]
hrt1a.0.3b <- pass1a.0.3b[[4]]
hrt1a.0.vars <- pass1a.0.vars[4,]
# save the data
save(hrt1a.0.nr1, hrt1a.0.1, hrt1a.0.nr2, hrt1a.0.3a, hrt1a.0.3b, hrt1a.0.vars,
  file = glue("{WD}/data/UM-rpneg/UM_rpneg_processed_hrt1a.0.RData"))

# kid
kid1a.0.nr1 <- pass1a.0.nr1[[5]]
kid1a.0.1 <- pass1a.0.1[[5]]
kid1a.0.nr2 <- pass1a.0.nr2[[5]]
kid1a.0.3a <- pass1a.0.3a[[5]]
kid1a.0.3b <- pass1a.0.3b[[5]]
kid1a.0.vars <- pass1a.0.vars[5,]
# save the data
save(kid1a.0.nr1, kid1a.0.1, kid1a.0.nr2, kid1a.0.3a, kid1a.0.3b, kid1a.0.vars,
  file = glue("{WD}/data/UM-rpneg/UM_rpneg_processed_kid1a.0.RData"))

# lun
lun1a.0.nr1 <- pass1a.0.nr1[[6]]
lun1a.0.1 <- pass1a.0.1[[6]]
lun1a.0.nr2 <- pass1a.0.nr2[[6]]
lun1a.0.3a <- pass1a.0.3a[[6]]
lun1a.0.3b <- pass1a.0.3b[[6]]
lun1a.0.vars <- pass1a.0.vars[6,]
# save the data
save(lun1a.0.nr1, lun1a.0.1, lun1a.0.nr2, lun1a.0.3a, lun1a.0.3b, lun1a.0.vars,
  file = glue("{WD}/data/UM-rpneg/UM_rpneg_processed_lun1a.0.RData"))

# liv
liv1a.0.nr1 <- pass1a.0.nr1[[7]]
liv1a.0.1 <- pass1a.0.1[[7]]
liv1a.0.nr2 <- pass1a.0.nr2[[7]]
liv1a.0.3a <- pass1a.0.3a[[7]]
liv1a.0.3b <- pass1a.0.3b[[7]]
liv1a.0.vars <- pass1a.0.vars[7,]
# save the data
save(liv1a.0.nr1, liv1a.0.1, liv1a.0.nr2, liv1a.0.3a, liv1a.0.3b, liv1a.0.vars,
  file = glue("{WD}/data/UM-rpneg/UM_rpneg_processed_liv1a.0.RData"))

# badi
badi1a.0.nr1 <- pass1a.0.nr1[[8]]
badi1a.0.1 <- pass1a.0.1[[8]]
badi1a.0.nr2 <- pass1a.0.nr2[[8]]
badi1a.0.3a <- pass1a.0.3a[[8]]
badi1a.0.3b <- pass1a.0.3b[[8]]
badi1a.0.vars <- pass1a.0.vars[8,]
# save the data
save(badi1a.0.nr1, badi1a.0.1, badi1a.0.nr2, badi1a.0.3a, badi1a.0.3b, badi1a.0.vars,
  file = glue("{WD}/data/UM-rpneg/UM_rpneg_processed_badi1a.0.RData"))

# wadi
wadi1a.0.nr1 <- pass1a.0.nr1[[9]]
wadi1a.0.1 <- pass1a.0.1[[9]]
wadi1a.0.nr2 <- pass1a.0.nr2[[9]]
wadi1a.0.3a <- pass1a.0.3a[[9]]
wadi1a.0.3b <- pass1a.0.3b[[9]]
wadi1a.0.vars <- pass1a.0.vars[9,]
# save the data
save(wadi1a.0.nr1, wadi1a.0.1, wadi1a.0.nr2, wadi1a.0.3a, wadi1a.0.3b, wadi1a.0.vars,
  file = glue("{WD}/data/UM-rpneg/UM_rpneg_processed_wadi1a.0.RData"))

# all tissues
save(pla1a.0.nr1, pla1a.0.1, pla1a.0.nr2, pla1a.0.3a, pla1a.0.3b, pla1a.0.vars,
     hip1a.0.nr1, hip1a.0.1, hip1a.0.nr2, hip1a.0.3a, hip1a.0.3b, hip1a.0.vars,
     gas1a.0.nr1, gas1a.0.1, gas1a.0.nr2, gas1a.0.3a, gas1a.0.3b, gas1a.0.vars,
     hrt1a.0.nr1, hrt1a.0.1, hrt1a.0.nr2, hrt1a.0.3a, hrt1a.0.3b, hrt1a.0.vars,
     kid1a.0.nr1, kid1a.0.1, kid1a.0.nr2, kid1a.0.3a, kid1a.0.3b, kid1a.0.vars,
     lun1a.0.nr1, lun1a.0.1, lun1a.0.nr2, lun1a.0.3a, lun1a.0.3b, lun1a.0.vars,
     liv1a.0.nr1, liv1a.0.1, liv1a.0.nr2, liv1a.0.3a, liv1a.0.3b, liv1a.0.vars,
     badi1a.0.nr1, badi1a.0.1, badi1a.0.nr2, badi1a.0.3a, badi1a.0.3b, badi1a.0.vars,
     wadi1a.0.nr1, wadi1a.0.1, wadi1a.0.nr2, wadi1a.0.3a, wadi1a.0.3b, wadi1a.0.vars,
  file = glue("{WD}/data/UM-rpneg/UM_rpneg_processed_pass1a.0.RData"))

8 Session Info

warnings()
session_info()
## ─ Session info ───────────────────────────────────────────────────────────────
##  setting  value                       
##  version  R version 4.0.2 (2020-06-22)
##  os       macOS  10.16                
##  system   x86_64, darwin17.0          
##  ui       X11                         
##  language (EN)                        
##  collate  en_US.UTF-8                 
##  ctype    en_US.UTF-8                 
##  tz       America/Toronto             
##  date     2021-12-15                  
## 
## ─ Packages ───────────────────────────────────────────────────────────────────
##  package      * version  date       lib source                                
##  assertthat     0.2.1    2019-03-21 [2] CRAN (R 4.0.2)                        
##  backports      1.2.1    2020-12-09 [2] CRAN (R 4.0.2)                        
##  broom          0.7.8    2021-06-24 [1] CRAN (R 4.0.2)                        
##  cachem         1.0.3    2021-02-04 [2] CRAN (R 4.0.2)                        
##  callr          3.7.0    2021-04-20 [1] CRAN (R 4.0.2)                        
##  cellranger     1.1.0    2016-07-27 [2] CRAN (R 4.0.2)                        
##  cli            3.0.1    2021-07-17 [1] CRAN (R 4.0.2)                        
##  coda           0.19-4   2020-09-30 [1] CRAN (R 4.0.2)                        
##  codetools      0.2-18   2020-11-04 [2] CRAN (R 4.0.2)                        
##  colorspace     2.0-2    2021-06-24 [1] CRAN (R 4.0.2)                        
##  crayon         1.4.1    2021-02-08 [2] CRAN (R 4.0.2)                        
##  curl           4.3.2    2021-06-23 [1] CRAN (R 4.0.2)                        
##  data.table     1.14.0   2021-02-21 [1] CRAN (R 4.0.2)                        
##  DBI            1.1.1    2021-01-15 [2] CRAN (R 4.0.2)                        
##  dbplyr         2.1.1    2021-04-06 [1] CRAN (R 4.0.2)                        
##  desc           1.3.0    2021-03-05 [1] CRAN (R 4.0.2)                        
##  devtools     * 2.4.2    2021-06-07 [1] CRAN (R 4.0.2)                        
##  digest         0.6.27   2020-10-24 [2] CRAN (R 4.0.2)                        
##  dplyr        * 1.0.7    2021-06-18 [1] CRAN (R 4.0.2)                        
##  ellipsis       0.3.2    2021-04-29 [1] CRAN (R 4.0.2)                        
##  evaluate       0.14     2019-05-28 [2] CRAN (R 4.0.1)                        
##  fansi          0.5.0    2021-05-25 [1] CRAN (R 4.0.2)                        
##  fastmap        1.1.0    2021-01-25 [2] CRAN (R 4.0.2)                        
##  forcats      * 0.5.1    2021-01-27 [2] CRAN (R 4.0.2)                        
##  fs             1.5.0    2020-07-31 [2] CRAN (R 4.0.2)                        
##  generics       0.1.0    2020-10-31 [2] CRAN (R 4.0.2)                        
##  ggfittext      0.9.1    2021-01-30 [2] CRAN (R 4.0.2)                        
##  ggplot2      * 3.3.5    2021-06-25 [1] CRAN (R 4.0.2)                        
##  glue         * 1.4.2    2020-08-27 [2] CRAN (R 4.0.2)                        
##  gridExtra      2.3      2017-09-09 [2] CRAN (R 4.0.2)                        
##  gtable         0.3.0    2019-03-25 [2] CRAN (R 4.0.2)                        
##  haven          2.3.1    2020-06-01 [2] CRAN (R 4.0.2)                        
##  highr          0.9      2021-04-16 [1] CRAN (R 4.0.2)                        
##  hms            1.1.0    2021-05-17 [1] CRAN (R 4.0.2)                        
##  htmltools      0.5.1.1  2021-01-22 [2] CRAN (R 4.0.2)                        
##  httr           1.4.2    2020-07-20 [2] CRAN (R 4.0.2)                        
##  impute       * 1.62.0   2020-04-27 [1] Bioconductor                          
##  inline         0.3.19   2021-05-31 [1] CRAN (R 4.0.2)                        
##  inspectdf      0.0.11   2021-04-02 [1] CRAN (R 4.0.2)                        
##  jsonlite       1.7.2    2020-12-09 [2] CRAN (R 4.0.2)                        
##  kableExtra   * 1.3.1    2020-10-22 [2] CRAN (R 4.0.2)                        
##  knitr          1.33     2021-04-24 [1] CRAN (R 4.0.2)                        
##  lattice        0.20-41  2020-04-02 [2] CRAN (R 4.0.2)                        
##  lifecycle      1.0.0    2021-02-15 [1] CRAN (R 4.0.2)                        
##  loo            2.4.1    2020-12-09 [1] CRAN (R 4.0.2)                        
##  lubridate      1.7.10   2021-02-26 [1] CRAN (R 4.0.2)                        
##  magrittr       2.0.1    2020-11-17 [2] CRAN (R 4.0.2)                        
##  MASS           7.3-53   2020-09-09 [2] CRAN (R 4.0.2)                        
##  matrixStats    0.60.0   2021-07-26 [1] CRAN (R 4.0.2)                        
##  memoise        2.0.0    2021-01-26 [2] CRAN (R 4.0.2)                        
##  modelr         0.1.8    2020-05-19 [2] CRAN (R 4.0.2)                        
##  MotrpacBicQC * 0.5.2    2021-06-25 [1] Github (MoTrPAC/MotrpacBicQC@43fb293) 
##  munsell        0.5.0    2018-06-12 [2] CRAN (R 4.0.2)                        
##  mvtnorm        1.1-2    2021-06-07 [1] CRAN (R 4.0.2)                        
##  naniar         0.6.1    2021-05-14 [1] CRAN (R 4.0.2)                        
##  pillar         1.6.2    2021-07-29 [1] CRAN (R 4.0.2)                        
##  pkgbuild       1.2.0    2020-12-15 [2] CRAN (R 4.0.2)                        
##  pkgconfig      2.0.3    2019-09-22 [2] CRAN (R 4.0.2)                        
##  pkgload        1.2.1    2021-04-06 [1] CRAN (R 4.0.2)                        
##  plyr           1.8.6    2020-03-03 [2] CRAN (R 4.0.2)                        
##  prettyunits    1.1.1    2020-01-24 [2] CRAN (R 4.0.2)                        
##  processx       3.5.2    2021-04-30 [1] CRAN (R 4.0.2)                        
##  progress       1.2.2    2019-05-16 [2] CRAN (R 4.0.2)                        
##  ps             1.6.0    2021-02-28 [1] CRAN (R 4.0.2)                        
##  purrr        * 0.3.4    2020-04-17 [2] CRAN (R 4.0.2)                        
##  R6             2.5.0    2020-10-28 [2] CRAN (R 4.0.2)                        
##  Rcpp           1.0.7    2021-07-07 [1] CRAN (R 4.0.2)                        
##  RcppParallel   5.1.4    2021-05-04 [1] CRAN (R 4.0.2)                        
##  readr        * 1.4.0    2020-10-05 [2] CRAN (R 4.0.2)                        
##  readxl         1.3.1    2019-03-13 [2] CRAN (R 4.0.2)                        
##  remotes        2.4.0    2021-06-02 [1] CRAN (R 4.0.2)                        
##  reprex         2.0.0    2021-04-02 [1] CRAN (R 4.0.2)                        
##  reshape2       1.4.4    2020-04-09 [2] CRAN (R 4.0.2)                        
##  rethinking   * 2.13     2021-08-13 [1] Github (rmcelreath/rethinking@2acf2fd)
##  rlang          0.4.11   2021-04-30 [1] CRAN (R 4.0.2)                        
##  rmarkdown      2.6      2020-12-14 [2] CRAN (R 4.0.2)                        
##  rprojroot      2.0.2    2020-11-15 [2] CRAN (R 4.0.2)                        
##  rstan        * 2.21.2   2020-07-27 [1] CRAN (R 4.0.2)                        
##  rstudioapi     0.13     2020-11-12 [2] CRAN (R 4.0.2)                        
##  rvest          1.0.0    2021-03-09 [1] CRAN (R 4.0.2)                        
##  scales         1.1.1    2020-05-11 [2] CRAN (R 4.0.2)                        
##  sessioninfo    1.1.1    2018-11-05 [2] CRAN (R 4.0.2)                        
##  shape          1.4.6    2021-05-19 [1] CRAN (R 4.0.2)                        
##  StanHeaders  * 2.21.0-7 2020-12-17 [1] CRAN (R 4.0.2)                        
##  stringi        1.7.3    2021-07-16 [1] CRAN (R 4.0.2)                        
##  stringr      * 1.4.0    2019-02-10 [2] CRAN (R 4.0.2)                        
##  testthat       3.0.4    2021-07-01 [1] CRAN (R 4.0.2)                        
##  tibble       * 3.1.3    2021-07-23 [1] CRAN (R 4.0.2)                        
##  tidyr        * 1.1.3    2021-03-03 [1] CRAN (R 4.0.2)                        
##  tidyselect     1.1.1    2021-04-30 [1] CRAN (R 4.0.2)                        
##  tidyverse    * 1.3.1    2021-04-15 [1] CRAN (R 4.0.2)                        
##  usethis      * 2.0.1    2021-02-10 [1] CRAN (R 4.0.2)                        
##  utf8           1.2.2    2021-07-24 [1] CRAN (R 4.0.2)                        
##  V8             3.4.2    2021-05-01 [1] CRAN (R 4.0.2)                        
##  vctrs          0.3.8    2021-04-29 [1] CRAN (R 4.0.2)                        
##  viridisLite    0.4.0    2021-04-13 [1] CRAN (R 4.0.2)                        
##  visdat         0.5.3    2019-02-15 [2] CRAN (R 4.0.2)                        
##  webshot        0.5.2    2019-11-22 [2] CRAN (R 4.0.2)                        
##  withr          2.4.2    2021-04-18 [1] CRAN (R 4.0.2)                        
##  xfun           0.24     2021-06-15 [1] CRAN (R 4.0.2)                        
##  xml2           1.3.2    2020-04-23 [2] CRAN (R 4.0.2)                        
##  yaml           2.2.1    2020-02-01 [2] CRAN (R 4.0.2)                        
## 
## [1] /Users/Alec/Library/R/4.0/library
## [2] /Library/Frameworks/R.framework/Versions/4.0/Resources/library